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In this paper the relationship between the density functional theory of freezing and phase field 
modeling is examined. More specifically a connection is made between the correlation functions that 
enter density functional theory and the free energy functionals used in phase field crystal modeling 
and standard models of binary alloys (i.e., regular solution model). To demonstrate the properties 
of the phase field crystal formalism a simple model of binary alloy crystallization is derived and 
shown to simultaneously model solidification, phase segregation, grain growth, elastic and plastic 
deformations in anisotropic systems with multiple crystal orientations on diffusive time scales. 
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I. INTRODUCTION 

The formalism for calculating equilibrium states was 
established many years ago by Gibbs, Boltzmann and 
others. While this formalism has proved remarkably suc- 
cessful there are many systems virhich never reach equilib- 
rium, mainly due to the existence of metastable or long 
lived transient states. This is most apparent in solid ma- 
terials. For example it is very unlikely that the reader is 
sitting a room containing any single crystals except items 
produced with considerable effort such as the silicon chips 
in computers. In fact the vast majority of naturally oc- 
curring or engineered materials are not in equilibrium 
and contain complex spatial structures on nanometer, 
micron or millimeter length scales. More importantly 
many material properties (electrical, optical, mechani- 
cal, etc.) are strongly influenced by the non-equilibrium 
structures that form during material processing. For ex- 
ample the yield strength of a polycrystal varies as the 
inverse square of the average grain size. 

The study of non-equilibrium microstructure forma- 
tion has seen considerable advances through the use 
^zjthe phase field approach. This methodology mod- 
els the dynamics of various continuum fields that col- 
lectively characterize microstructure in phase transfor- 
mations. For example, phase field or continuum models 
have been used to simulate spinodal decomposition (Ij], 
order-disorder transition kinetics ordering of block- 
copolymer melts Q , solidification of pure and binary sys- 
tems P,|EIS0|8l and many other systems. In these phe- 
nomena the evolution of the appropriate field(s) (e.g., so- 
lute concentration in spinodal decomposition) is assumed 
to be dissipative and driven by minimizing a phenomeno- 
logical free energy functional J^J . 

Advances in the phase field modeling of solidifica- 
tion phenomena have followed a progression of inno- 
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vations, beginning with the development of free ener- 
gies that capture the thermodynamics of pure materi- 
als 0, 0, 01 and alloys Hi . Several modification were 
then proposed llOl lllj| to simplify numerical simu- 
lations and improve computational efficiency. Perhaps 
the most important innovation was the development of 
matched asymptotic analysis techniques that directly 
connect phase field model parameters with the classi- 
cal Stefan (or sharp-interface) models for pure materials 
or alloys 0, 0, 0|. These techniques were com- 
plimented by new adaptive mesh refinement algorithms 
[l6Lll7l |. whose improved efficiency significantly increased 
the length scales accessible by numerical simulations, 
thus enabling the first experimentally relevant simula- 
tions of complex dendritic structures and their interac- 
tions in organic and metallic alloys 18, JJ,,^, ^1, 22]. 

A weakness of the traditional phase field methodology 
is that it is usually formulated in terms of fields that are 
spatially uniform in equilibrium. This eliminates many 
physical features that arise due to the periodic nature of 
crystalline phases, including elastic and plastic deforma- 
tion, anisotropy and multiple orientations. To circum- 
vent this problem traditional phase field models have 
been augmented by the addition of one or more auxil- 
iary fields used to describe the density of dislocations 
Ea El i3 continuum stress and strain fields [H E^ 
and orientation fields EEESE3- These approaches have 
proven quite useful in various ap plications such as poly- 
crystalline solidification 0, Effl, 29, El E3- Nev- 
ertheless it has proven quite challenging to incorporate 
elasto-plasticity, diffusive phase transformation kinetics 
and anisotropic surface energy effects into a single, ther- 
modynamically consistent model. 

Very recently a new extension to phase field model- 
ing has emerged known as the phase field crystal method 
(PFC) iSS, J54, 35] . This methodology describes the evo- 
lution of the atomic density of a system according to 
dissipative dynamics driven by free energy minimization. 
In the PFC approach the free energy functional of a solid 
phase is minimized when the density field is periodic. As 
discussed in prior pubhcations [s^, 0, E3 the periodic 
nature of the density field naturally gives rise to elastic 
effects, multiple crystal orientations and the nucleation 
and motion of dislocations. While these physical fea- 
tures are included in other atomistic approaches (such as 
molecular dynamics) a significant advantage of the PFC 
method is that, by construction, it is restricted to oper- 
ate on diffusive time scales not on the prohibitively small 
time scales associated with atomic lattice vibrations. A 
similar approach has also been recently proposed by Jin 
and Khachaturyan |3Q] . In the case of pure materials the 
PFC approach has been shown [s^, E3 to model many 
phenomena dominated by atomic scale elastic and plas- 
tic deformation effects. These include grain boundary 
interactions, epitaxial growth and the yield strength of 
nano-crystals. 

The original PFC model is among the simplest mathe- 
matical descriptions that can self-consistently combine 



the physics of atomic-scale elasto-plasticity with the 
diffusive dynamics of phase transformations and mi- 
crostructure formation. Nevertheless, analogously to tra- 
ditional phase field modeling of solidification, further 
work is required to fully exploit the methodology. More 
specifically it is important to be able to generalize the 
method to more complex situations (binary alloys, faster 
dynamics, different crystal structures, etc.), to develop 
more efficient numerical techniques and to make a di- 
rect connection of the parameters of the model to ex- 
perimental systems. Several innovations toward this goal 
have already been developed. Goldenfeld et al. 1401] 
have recently derived amplitude equations for the PFC 
model which are amenable to adaptive mesh refinement 
schemes. This work has the potential to enable simu- 
lations of mesoscopic phenomena (/im mm) that are 
resolved down to the atomic scale and incorporate all the 
physics discussed above. Another recent advance is the 
inclusion of higher order time derivatives in the dynamics 
to simulate "instantaneous" elastic relaxation |41|. This 
extension is important for modeling complex stress prop- 
agation and externally imposed strains. Very recently, 
Wu et al. [sl] fitted the PFC parameters to experimen- 
tal data in iron and were able to to show that the PFC 
model gives an accurate description of the anisotropy of 
the surface tension. In addition to this work Wu and 
Karma have also developed a simple and elegant scheme 
to extend the method to other crystal symmetries (i.e., 
FCC in three dimensions). 

The purpose of this paper is to link the formalism of 
density function theory (DFT) of freezing, as formulated 
by Ramakrishnan and Yussouff (and also reviewed by 
many other authors, such as Singh |5^) with the phase 
field crystal (PFC) method and to exploit this connection 
to develop a phase field crystal model for binary alloys. 
The organization of the paper and a summary of the 
remaining sections is as follows. 

In Section IIA the density functional theory of freez- 
ing of pure and binary systems is briefly outlined. In this 
approach the free energy functional is written in terms 
of the time averaged atomic density field p {pA and pB 
in binary systems) and expanded around a liquid refer- 
ence state existing along the liquid/solid coexistence line. 
Formally the expansion contains the n-point correlation 
functions of the liquid state. In this work the series ex- 
pansion of the free energy is truncated at the 2-point 
correlation function, C{fi,f2)- 

Within this framework it is shown in Section IIIA that 
the PFC model for a pure material can be recovered from 
DFT if C'^fi, r2) is parameterized by three constants re- 
lated to the liquid and solid state compressibilities and 
the lattice constant. The parameters of the PFC model 
can thus be directly related to the physical constants that 
enter the DFT of freezing and the PFC model can be 
viewed as a simplified form of DFT. In Section IIIB a 
binary system is considered. Similar to the case of pure 
materials the free energy expansion of a binary alloy will 
be truncated at the 2-point correlation functions which 
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are then characterized by three parameters. At this level 
of simplification it is shown that the "regular" solution 
model used in materials physics for alloys can be obtained 
directly from density functional theory. It is shown that 
the phenomenological nearest neighbour bond energies 
that enter the "regular" solution model are equal to 
the compressibilities that enter density functional theory. 
This section also provides insight into the concentration 
dependence of various properties of the crystalline phase 
of a binary alloy such as the lattice constant, effective 
mobilities and elastic constants. 

In Section IIIC a simplified version of the binary alloy 
free energy is derived. This is done in order to provide 
a mathematically simpler model that can more transpar- 
ently illustrate the use of the PFC formalism in simulta- 
neously modeling diverse processes such as solidification, 
grain growth, defect nucleation, phase segregation and 
elastic and plastic deformation. This section also shows 
that the free energy of the simplified alloy PFC model 
reproduces two common phase diagrams associated with 
typical binary alloys in materials science. 

In Section IV dynamical equations of motion that gov- 
ern the evolution of the solute concentration and density 
field of the binary alloy are derived. Finally in Section 
V the simplified binary alloy model is used to simulate 
several important applications involving the interplay of 
phase transformation kinetics and elastic and plastic ef- 
fects. This includes solidification, epitaxial growth and 
spinodal decomposition. Some of the more tedious cal- 
culations are relegated to the Appendices. 



II. DENSITY FUNCTIONAL THEORY OF 
FREEZING 

In this section free energy functionals of pure and bi- 
nary systems as derived from the density functional the- 
ory of freezing are presented. A derivation of the func- 
tional for a pure materials is outlined in the Appendix 
IbI For a more rigorous treatment the reader is referred 
to the work of Ramakrishnan and Yussouff 52] and nu- 
merous other very closely related review articles by Singh 
Eal, Evans Ej and references therein. 



A. Single Component System 

In density functional theory the emergence of an or- 
dered phase during solidification can be viewed as a tran- 
sition to a phase in which the atomic number density, 
/9(r), is highly non- homogenous and possesses the spatial 
symmetries of the crystal 53J. This approach implic- 
itly integrates out phonon modes in favour of a statisti- 
cal view of the ordered phase that changes on diffusive 
time scales. The free energy functional of a system is ex- 
pressed in terms of p and constitutes the starting point 
of the PFC model. 



In this work the free energy functional, denoted T[p\, 
is expanded functionally about a density, p = pi, corre- 
sponding to a liquid state lying on the liquidus line of the 
solid-liquid coexistence phase diagram of a pure material 
as shown in Fig. (^). The expansion is performed in 
powers oi Sp = p ^ pi. 

As shown by other others and outlined in Appendix IbI 
the free energy density can be written as 



knT 



dx 



p{r) In ( ] - Sp{r) 



Pi 



~2 I dridr25p{ri)C2{ri,r2)Sp{r2) + 
~\ I dfidr2dr3^p{ri)C3{ri,r2,r3)6p{f2)Sp{'r3) 
+ ••• (1) 

where J-c is the free energy corresponding to the density 
p{r) minus that at the constant density pi. The func- 
tion, C2 is the two point direct correlation function of an 
isotropic fluid. It is usually denoted dj = (72(^12) and 
satisfies ri2 = [ri — r2|- The function C3 is the three 
point correlation function, etc.. Formally the correlation 
functions are defined by 



C2(ri,r2) 

C3(ri,r2,f3) 



Sp{r) 



dp{fi)Sp{r2) 
53$ 



Sp{r3)Sp{ri)Sp{r2) 



(2) 



where $[p] represents the total potential energy of inter- 
actions between the particles in the material. The proof 
that $ is a functional of p is shown rigorously in Evans 



B. Multi-Component System 

For an alloy involving one or more components the 
free energy functional of a pure material in Eq. is 
extended to the form 



E 



df 



Pt {r} In {pt (r) / p\) - 5p^ {r) 



- j dri dr2 Spi (fi ) Cy {n , r2 ) (r2 ) H 

(3) 

where the sums are over the elements in the alloy, Spi = 
Pi — Pi, Pg is the value of the number density of compo- 
nent i on the liquid side of the liquid/solid coexistence 
line. The function Cy is the two point direct correlation 
function of between components i and j in an isotropic 
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given in Eq. Q can be truncated at C2, i.e., 



C(k) 




FIG. 1: (a) Sample phase diagram. In this figure the shaded 
area corresponds to a coexistence region. In the calculations 
presented in this paper the correlation functions are taken 
from points along the liquidus line at density p^. (b) In this 
figure a 'typical' liquid state two-point direct correlation func- 
tion is sketched. The dashed line represents the approxima- 
tion used in most of the manuscript. 



fluid. As in the case of a pure materials it will be assumed 
that Cij = Cijir 12), where ri2 = 1^1—^21. The next term 
in the expansion of Eq. Q contains the three point cor- 
relation, the next after that, the four point, etc.. In this 
paper only two point correlations will be considered, but 
it must be stressed that these higher order correlations 
maybe crucial for some systems, such as Si. Before con- 
sidering the properties of the binary alloy free energy in 
detail it is instructive to first study the properties of a 
pure system and show the connection between this for- 
malism and the original phase field crystal model. 



III. ANALYSIS OF FREE ENERGY 
FUNCTIONALS 



Pure Materials 



In this section the free energy functional of a single 
component alloy is considered in the limit that the series 



T/ksT = j dr [p\n{p/ pi) - 5p\ 

- (1/2) / dfdf'5pC{r,r')5p' (4) 



where, for convenience, the subscript '2' has been 
dropped from the two-point correlation function as has 
the subscript c from Tc- To understand the basic fea- 
tures of this free energy functional it is useful to expand 
T in the dimensionless deviation of the density p from 
its average, p, using the re-scaled density 



= {p~ p)Ip- 

Expanding J- in powers of n gives. 



pksT 



= / dr 



1 



- pC n n 

"2— "-T + T2 



(5) 



, (6) 



where IS.J- = T — To and To is the the free energy func- 
tional at constant density (i.e., p — p). For simplicity C 
is an operator defined such that nCn = J dr 'n(r)C(|r — 
r'l)n(r'). Terms that are linearly proportional to n in 
the above integral arc identically zero by definition. 

To gain insight into the properties of the free energy 
functional in Eq. © it is useful to expand the two point 
correlation function in a fourier series, i.e., 

C^Cq + C2k^ + Cik^ + • • • . (7) 

(in real space this corresponds to C = {Cq — C^SI^ + 
6*4 V* _ . . . ),5(j?_ 7^)^ where the gradients are with re- 
spect to r'). The function C is sketched for a typical 
liquid in Fig. ^p)- In what follows only terms up to 
A:^ will be retained. In this manner the properties of the 
material are parameterized by the three variable, Co, Ci 
and C4. To fit the first peak in (7, Cq, C2 and must be 
negative, positive and negative, respectively. These vari- 
ables are related to three basic properties of the material, 
the liquid phase isothermal compressibility (^ (1 — pCo)), 
the bulk modulus of the crystal (~ pG\l\G/S) and lat- 
tice constant (~ {C2I\Ca\)^I'^). In other words the /c = 
term is related to the liquid phase isothermal compress- 
ibility, the height of the first peak (Cm in Fig. (^)) 
is related to the bulk modulus of the crystalline phase 
and the position of the first peak determines the lattice 
constant. 

It is important to note that at this level of simplifica- 
tion the material is only defined by three quantities which 
may not be enough to fully parameterize any given ma- 
terial. For example this simple three parameter model 
always predicts triangular symmetry in two dimensions 
and BCC symmetry in three dimensions. Other crystal 
symmetries can be obtained by using more complicated 
2-point correlation functions |5l| or by including higher 
order correlation functions. 
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In two dimensions T is minimized by a triangular lat- 
tice that can be represented to lowest order by a one- 
mode approximation as 



Substituting Eq. ^ into Eq. © gives, 



(8) 



AF = -(l-pCo)^ 



2 1 ^3 ^ i5 ^4 _ . . . 



32 



512 



A' 



(9) 



where AF = AT j {pkBT S) and S is the area of a unit 
cell. Minimizing AF with respect to q gives the equilib- 
rium wavevector, q^q, which is 



3(72/(8|C4|). 



(10) 



or in terms of the equilibrium lattice constant aeq ~ 
2tt /qeq. When q — q^q, AT becomes 



AF = ^ABA^--A'' + — 
16 32 512 



(11) 



where AB = B^ - B^ = 1 - pCo and = 
P (C'2)^/(4|C'4|). The parameter B^ is the dimensionless 
bulk modulus of the liquid state (i.e., B^ = n/ {pkBT), 
where n is the bulk modulus of a liquid). The parameter 
Bs is proportional to the bulk modulus in the crystalline 
phase. 

Equation pi|l indicates that the liquid state is linearly 
unstable to the formation of the crystalline phase when 
AB < 0. This instability arises from a competition be- 
tween the elastic energy stored in the liquid and crys- 
talline phases. It is interesting to note that AB can also 
be written, 



AB = (p, - p)/ps 



(12) 



where ps = 1 /Cm and C„i is the height of the first peak 
of C as shown in Fig. Written in this form p^ can be 
thought of as defining the effective spinodal density, i.e., 
the average density at which the liquid becomes linearly 
unstable to crystallization. 

Unfortunately it is difficult to obtain the equilibrium 
state (i.e., by solving dAF/dA — 0) without truncating 
the infinite series in Eq. lll|l . If only terms to order A* 
are retained an analytic approximation can be obtained 
for the amplitude (Amm) that minimizes T. In this ap- 
proximation the solution is 

(l + v/20p/p, -19)/5. (13) 

Thus solutions for a crystalline state exist when p > 
19/20 p,. 

It is also straightforward to calculate the change in en- 
ergy of the crystalline state upon deformation (i.e., bulk. 



shear or deviatoric). Details of similar calculations are 
given in previous pubhcations [s^, 13- The result of 
these two-dimensional calculations gives the dimension- 
less bulk modulus. Be, of the crystalline phase, i.e.. 



(14) 



In other words the parameter controls the bulk mod- 
ulus of the crystalline phase. 

These calculations can be easily extended to three di- 
mensions. As discussed previously this particular ap- 
proximation for C leads to a BCC crystal in three di- 
mensions which can be represented in a one mode ap- 
proximation as, 

n = ^( cos (gx) cos ((/y) + cos ((/a;) cos (qz) 

+ cos(g2/)cos((7z)). (15) 

Substituting this functional form into the free energy and 
minimizing with respect to q gives. 



<=V^2/|C4|, 

and the free energy functional at this q is, 

AF^'' = -ABA' - -A^ + —A^ + 
8 8 256 



(16) 



(17) 



Thus in this instance the 'spinodal' occurs at the same 
density as in the two dimensional case. If the series is 
truncated at the amplitude that minimizes the free 
energy is then 



A 



3d 



4 (1 + ^15f5/f5s - 14) / 



15. 



(18) 



Thus in this approximations crystalline (BCC) solutions 
only exist if p > lAps- In addition the elastic constants 
can also be calculated in 3D in the usual manner. For ex- 
ample the dimensionless bulk modulus of the crystalline 
state is given by; 



3d 



(19) 



This calculation gives the basic functional dependence 
of the (dimensionless) bulk modulus on B'^ and the am- 
plitude. For a more accurate calculation higher order 
fourier components and more terms in the powers series 
in n should be retained. 

Finally it is useful to consider fixing the density and 
varying the temperature. If the liquidus and solidus lines 
are roughly linear then, ps can be approximated by a 
linear function of temperature. In the sample phase di- 
agram shown in Fig. (^) the liquidus and solidus lines 
are roughly parallel and it is likely that the spinodal is 
also roughly parallel to these lines. In this case AB can 
be written, 



AB = QfAT 



(20) 



where AT = {T — Tg) /Ts, Tg is the spinodal temperature 
and a = [Ts/ Ps){dps/ dTg), evaluated dX p = p. 
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B. Binary Alloys 

For a binary alloy made up of and 'B' atoms the 
free energy functional can be written to lowest order in 
terms of the direct correlation functions as, 



= / dr 

1 
2 



PA In 



- Sp/ 



PB In ( ^ j - 6pi 



dfi dr2 5pA [ri ) Caa {n , f2 ) SpA {r2 ) 
+SpBiri) CBB{ri,r2) 5pB{r2) 
+2 5pA{ri) Cab (7^1,7^2) 5pB{r2) 



(21) 



where 5pA = Pa — pf and 6pB = Pb — pf ■ It is assumed 
here that all two point correlation functions are isotropic, 
i.e., Ctj{fi,f2) = Cij{\fi -r2\)- 

In order to make a connection between the alloy free 
energy and standard phase field models it is useful to 
define the total number density, p = pa + Pb and a local 
concentration field c = pa/ P- In terms of these fields 
the atomic densities can be written, pA = cp and pb = 
p{l — c). Furthermore it is useful to define p = pi + Sp 
where pi = pf + pf and Sc = 1/2 — c. Substituting these 
definitions into Eq. (|^ gives, 

= f df p\n{p/pi) - 6p ^ 



1 



--Sp{cCaa 



[1-c)Cbb]Sp 
+/9{cln (c) + (1 - c) ln(l - c) } 
+pc/\C{l-c)p + p5c + Fo 



where 



AC 



{Caa 



AA 



Cbb)/2 — Cab, 

- Cbb){p + Pi) + P In 



pI 



(22) 



(23) 



and 




Caa 



{(pfr + '^ipi+p)) 



Cbb 



(ipfr + ^ipi+p)) 



(24) 



To illustrate the properties of the model in Eq. 1)22(1 
it useful to consider two limiting cases, a liquid phase at 
constant density and a crystalline phase at constant con- 
centration. These calculations are presented in following 
two sub-sections. 



1. Liquid phase properties 

In the liquid phase p is constant on average and in the 
mean field limit can be replaced hy p = p. To simplicity 



calculations, the case p = p ~ pi [or 5p ~ 0) will now 
be considered. As in the previous section it is useful to 
expand the direct correlation functions in Fourier space, 



I.e. 



Ci j — Cq -\- C2 k c ^ k -\- ' ' ■ 



(25) 



where the subscript i and j refer to a particular element. 
Substituting the real-space counterpart of the Fourier ex- 
pansion for Cij (to order k'^) into Eq. (|22|l gives, 



pksT 



dr 
pACo 



cln (c) + (1 - c) ln(l - c) 



c(l 



c)+/fc+^|Vcp 



,(26) 



where J^c is the total free energy minus a constant that 
that depends only on p, and p'^. 



— / dBB 



Bt^)+-p\n{pf/pf) 



ACn = C 



_ AAA 



AB 



(27) 



(28) 



and By = 1 — pCq is the dimensionless bulk compress- 
ibility. Equation H2t)|l is the regular solution model of a 
binary alloy. 

The coefficient of c(l — c) in Eq. (|26|) is given by 



pACo = 2Bf'' - Bf^ ~ Bf"" . 



(29) 



This result shows that in the liquid state the 'interaction' 
energies that enter regular solution free energies are sim- 
ply the compressibilities (or the elastic energy) associated 
with the atomic species. The 7^ term is also quite inter- 
esting as it is responsible for asymmetries in the phase 
diagram. Thus Eq. (|27|) implies that asymmetries can 
arise from either different compressibilities or different 
densities. 

Expanding Eq. (|26|l around c = 1/2 gives. 



ATc 
pksT 



dr 



2 



4 



(30) 



where, AJ^c = ~ pksT f dr{pACo/8 - ln(2)), u = 
16/3, / = (4 - pACo) and K = pAC2. The parameter 
is related only to the fc = part of the two-point 
correlation function and can be written. 



AA 



B 



BB 



2Br) 



(31) 



This result implies that the instability to phase segre- 
gation in the fluid is a competition between entropy (4) 
and the elastic energy of a mixed fluid {2Bf^) with the 
elastic energy associated with a phase separated fluid 
{Bf^ + Bf^). Replacing the dimensionless bulk moduli 
with the dimensional version (i.e., B — n/kBT), gives 



the critical point (i.e. 



0) as 



= {2k AB - KBB - KAA)/{4:kB)- 



(32) 
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where k is the dimensional bulk modulus. 

The properties of the crystalline phase are more com- 
plicated but at the simplest level the only real difference 
is that the elastic energy associated with the crystalline 
state must be incorporated. This is discussed in the next 
section. 



2. Crystalline phase properties 

To illustrate the properties of the crystalline state, the 
case in which the concentration field is a constant is con- 
sidered. In this limit the free energy functional given in 
Eq. (|22l) can be written in the form 



= / dr 



pin ( -^1 -5p-]-5pC5p + G 



(33) 



where G is a function of the concentration c and pi and 
couples only linearly to bp. The operator C can be writ- 
ten as 



C 



^CAA + {l-cfCBB + M^~c)CAB- (34) 



Thus in the limit that the concentration is constant this 
free energy functional is that of a pure material with 
an effective two point correlation function that is an av- 
erage over the AA, BB and AB interactions. In this 
limit the calculations presented in section IIII Al can be 
repeated using the same approximations (i.e., expand- 
ing p around p^, expanding C to and using a one 
mode approximation for 5p) to obtain predictions for 
the concentration dependence of various quantities. For 
example the concentration dependence of the equilib- 
rium wavevector (or lattice constant, Eq. (|10|) ') and 
bulk modulus Eq. (|14|) can be obtained by substituting 
(7„ = c^C^^ + {l -cfC^^ + 2c(l - c)C^B. 

As a more specific example the equilibrium lattice con- 
stant can be expanded around c = 1/2 to obtain in two 
or three dimensions, 



aeq{5c) = aeq(O) {I + T] 5c 



(35) 



where 5c — c — 1/2 and r/ is the solute expansion coeffi- 
cient given by, 



where 



77 = {5Ca - 5C2)/2 



(36) 



(37) 



and (7° ^ Cn{5c = 0)^ - ((7,f ^ + C^b + 2C'^«)/4. 

This line of reasoning can also be used to understand 
the influence of alloy concentration on crystallization. 
Specifically, for the case of an alloy, the terms in Eq. 

(with A replaced with Amin) become functions of 
concentration, since AS and Amin are concentration de- 
pendent. Here, AB can be expanded around c — 1/2, 



I.e., 



AB{Sc) = ABq + ABi 5c + AB2 Sc^ 



(38) 



where ABq 



Sg, ABi = B{- Bf and AB2 



B2 ~ B2 arc determined in the appendix. This would 
imply that in the crystalline phase the free energy has a 
term of the form, r'^((5c)^, where 

3AB2 ,2 



pACo(l-H3AVj8) 



3/8i?|AV„, (39) 



in two dimensions (in three dimensions the 3/8 factor is 
replaced with 3/4). This result indicates that crystalliza- 
tion (i.e., a non-zero A^m) favours phase segregation, as- 
suming kaa + kbb < '^kab- For example, when i3| = 0, 
the critical temperature increases and can be written, 

r5 = TS(l + 3AL„/8), (40) 

or ^TIj{1 + 3v4,2„j„/4) in three dimensions. 

C. Simple Binary Alloy Model 

In this section a simple binary alloy model is pro- 
posed based on a simplification of the free energy in 
Eq. H22|l . The goal of this section is to develop a math- 
ematically simple model that can be used to simultane- 
ously model grain growth, solidification, phase segrega- 
tion in the presence of elastic and plastic deformation. To 
simplify calculations it is convenient to first introduce the 
following dimensionless fields. 



nA 
hb 



[PA 
{PB 



pa)Ip 
pb)Ip- 



(41) 



Also, it is convenient to expand in the following two fields. 



5N 



nA + riB 
[riB -nA) + 



PB - PA 



(42) 



The following calculations will use the field 5N instead 
of 5c. Expanding Eq. around 5N = and n = 

gives a free energy of the form 



pUbT 



J dr{^ [B, + B, {2R'W' + R^W')] n 

^ - ^n^ + ^f5N+-5N^ + -5N'' 
4 2 4 



(43) 



Details of this free energy and explicit expressions for 
B^ , B^ and R are given in Appendix |0 For simplicity 
the calculations presented in this section are for a two 
dimensional system. 

The transition from liquid to solid is intimately related 
to AB ~ B^ ~ B" as was the case for the pure material 
and can be written in terms of a temperature difference, 
i.e., Eq. (|20|l . In addition some of the polynomial terms 
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in n and 6N have been multiplied by variable coefficients 
even though they can be derived exactly as shown in 
Appendix [0 For example the parameter v — 1/3 re- 
covers the exact form of the term. This flexibility in 
the choice of coefficients was done to be able to match 
the parameter of the free energy with experimental ma- 
terials parameters. As an example Wu and Karma [sH 
showed that adjusting the parameter v can be used to 
match the amplitude of fluctuations obtained in molecu- 
lar dynamics simulations. With this fit they are able to 
accurately predict the anisotropy of the surface energy of 
a liquid/crystal interface in iron. 

To facilitate the calculation of the lowest order phase 
diagram corresponding to Eq. (|43|l it is convenient to as- 
sume the concentration field 5N varies significantly over 
length scales much larger than the atomic number density 
field n. As a result, the density field can be integrated 
out of the free energy functional. In this instance the 
one-mode approximation for n, i.e., 



/ 1 f2qy\ qy 
n = A \ - cos —pr — cos(qx) cos — ^ — 



(44) 



wiU be used. Substituting Eq. into Eq. lO and 

minimizing with respect to q and A (recalling that SN is 
assumed constant over the scale that n varies) gives 



V3/(2i?) 



(45) 



and 



AminiSN = 0) (which is thus a function 
of ABq). This simple form can be used to calculate 
the solid/solid coexistence concentrations at low temper- 
atures according to 



(50) 



The critical temperature, ABq is determined by setting 
SNcoex ~ and solving for ABq, which gives. 



AB^ = 15wv - 2tJ^6B'^w / (65^ 



(51) 



2. Solid-liquid coexistence 

To obtain the liquid/solid coexistence lines the free 
energy of the liquid state must be compared to that of 
the solid. The mean field free energy of the liquid state 
is obtained by setting n = which gives. 



Itq 2 4 



(52) 



To obtain the solid-liquid coexistence lines it is useful 
to expand the free energy of the liquid and solid states 
around the value of 6N at which the liquid and solid 
states have the same free energy, i.e., when Fsoi = Fnq. 
This occurs when, 



and 



4 
15 



-(t + y/t^ - 15wAb) . 



(46) 



The free energy that is minimized with respect to ampli- 
tude and lattice constant is then, 



> 

45w 
'512' 



16 



-ABAi 



t 

16' 



-At 



(47) 



For mathematical simplicity all further calculations 
will be limited to the approximations B^ = Bq + B2{SN)'^ 
and B'^ = Bq. In this limit analytic expression can be 
obtained for a number of quantities and the free energy 
functional is still general enough to produce for example 
a eutectic phase diagram. 



1. Solid-solid coexistence 

Expanding Fgoi around SN — gives, 
Fsoi=FsoiiO)+a6Ny2 + b6N^/4 
where 

a EE w + 3Bi{A°^,jy8 



(48) 



\2 AO 

min 



(15z;<,„ - At) 



(49) 



5N,, 



± 



[ABl^ 



ABo) /Bi, 



(53) 



where AB^ = B'q - B'q and AB'q^ = 8tV(135w) is the 
lowest value of AB\^ at which a liquid can coexist with 
a solid. To complete the calculations, Fgoi and Fuq are 
expanded around 5Nis to order {SNis)-, i.e.. 



Fug = TiSNis) + au,{SN ~ SNis) +bug{6N ^ SNis)' 
Fsoi = J^iNis) + a,oi{5N-SNis)+bsoi{5N-5Nis)^ 



(54) 



(note by definition Fsoi{6Nis) = Fiiq{6Nis)), where 



aiiq 
O-sol 
hiq 



{w + u6Nl)5Nu 

auq + m'^BiSNis/iQlbv^) 

[w + 3uSN?,)/2 



6 

5v 



B' 



-ABq 



AB'q' 



(55) 



It is now straightforward to calculate the liquid/solid co- 
existence lines from Eq. I|54() in terms of the parameters 
given above. The liquid/solidus lines are; 



SNi 



iq 



SNsol = 




.(56) 
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for SNiiq > 0, SNsoi > and similar results for SNug < 0, 
SNsoi < 0, since F is a function of SN"^ in this example. 
The calculations in this section and the previous section 
are reasonably accurate when ABq^ > ABg, however in 
the opposite limit a eutectic phase diagram forms and 
the accuracy of the calculations decreases. This case will 
be discussed below. 



3. Linear clastic constants 

As shown in previous publications [s^, 0| , the elastic 
constants can be calculated analytically in a one mode 
approximation by considering changes in _F as a function 
of strain. For the binary model similar calculations can 
be made in a constant SN approximation and give, 

Cii/3 = Ci2 = Cm = ^B'{A^^nf. (57) 

(this calculation can be done for arbitrary dependence of 
on SN). As expected the elastic constants are directly 
proportional to the amplitude of the density fluctuations. 
This implies that the elastic constants decrease as the 
liquid solid transition is approached from the solid phase. 
This result implies both a temperature and concentration 
dependence through the dependence of Amin on AB. In 
addition to this dependence (which might be considered 
as a 'thermodynamic' dependence) the magnitude of the 
elastic constants can be altered by the constant B'' . 

4- Calculation of phase diagrams 

To examine the validity of some of the approxima- 
tions made in the previous section, numerical simulations 
were conducted to determine the properties of the solid 
and liquid equilibrium states. The simulations were per- 
formed over a range of SN values, three values of ABq 
and two values of w. The specific values of the constants 
that enter the model are given in the figure captions. 

In Fig. ^ analytic predictions and numerical solu- 
tions are shown for the free energy, F, the lattice constant 
R and bulk modulus at three values of ABq, using the 
free energy functional given in Eq. I|43|l . These figures 
indicate that the approximate solutions are quite accu- 
rate except for some deviations in the SN dependence 
of the bulk modulus. For the specific set of parameters 
used for these comparisons, the one-mode constant con- 
centration approximation predicts no asymmetry in any 
quantity. However, it is clear from the numerical solu- 
tions that some symmetry does exist. Figure Q shows 
the phase diagram corresponding to the parameters used 
in Fig. 

The same calculations presented in Figs. Q and © 
were repeated using w — -0.04 for which AB^ > AB^". 
In this case the agreement between the numerical results 
and the one-mode constant concentration calculations for 
the free energy, lattice constant and bulk modulus are 




(5N 



FIG. 2: Free energy (a), lattice constant (b) and bulk modulus 
(c) are plotted as a function of SN for three different values of 
ABq. In figure (a) the lines from top to bottom (and bottom 
to top in (c)) correspond to ABq = 0.07, 0.02 and -0.03 
respectively computed in the one-mode approximations. The 
solid, starred and open points correspond to ABo = 0.07, 
0.02 and —0.03 and were calculated directly by numerically 
by minimizing the free energy functional given in Eq. (1431 . 
Other parameters used correspond Bo = 1.00, Bi — 0, B2 = 
-1.80, t = 0.60, V = 1.00, w = 0. 088, u = 4.00, L = 4.00, 
Ro = 1.00 and Ri = 1.20 (see Eg. 1^121 for definitions of Ro 
and Ri). 



similar to the w = 0.088 case as shown in Fig. (0}. The 
phase diagram corresponding to the parameters used in 
Fig. Q is shown in Fig. jsj. In this case, the analytic cal- 
culations (Eqns. (|56|l and H50|) ') for the coexistence lines 
breakdown at the eutectic point. Higher order terms in 
the in SN are needed to accurately predict the coexis- 
tence lines. 



IV. DYNAMICS 

To simulate microstructure formation in binary alloys, 
dynamical equations of motions for the field SN and n 
need to be developed. The starting point is the full free 
energy in Eq. lj2U, written in terms of in terms of pA 
and pb, i.e., !F{pA, Pb)- The dynamics of pA and pB 
is assumed to be dissipative and driven by free energy 
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FIG. 3: Phase diagram of ABo Vs. 5N for the parameters 
corresponding to the parameters of Fig. The sohd hne 

is a numerical solution of the one mode approximation and 
the dashed lines are from Eq. 15UII for the lower solid/solid 
coexistence lines and Eq. 15611 for the upper liquid/solid co- 
existence lines. The solid points are from numerical solutions 
for the minimum free energy functional given in Eq. 14311 . 



FIG. 5: Phase diagram of ABo Vs. 5N for the same parame- 
ters as those used to generate Fig. ^ , with the exception that 
w = —0.04. The dotted lines below the eutectic temperature, 
correspond to metastable states. 



AB^ 0.028 



minimization, i.e., 




FIG. 4: This figure is identical to Fig. Q with the exception 
that w = -0.04. 



dpA 
dt 

dpB 
dt 



Ma (pa.Pb) V 
Ms (pa, PS )V 



5J^ 

SpA 

SpbJ 



(58) 
(59) 



where Ma and AIb are the mobihties of each density field. 
In general these mobihties depend on the density of each 
species. The free energy T{pa,Pb) can equivalently be 
defined in terms of n and 5N . This allows the previous 
equations to be re-written as 

,%i.(M,v',v^M.v)f|£-^) m 



dt 
-dpB 



{MbV^ + V • MbV) ( — 



my) ■ 



Adding Eq. H60(l and Eq. 16111 and collecting terms gives 
the following equation for n, 



dn 
'dt 



{A/iV^ + VMi • V} 



5T 



{MaV^ + VM2 • V} 



d{6N) ' 



(62) 



where Mi = {Ma + Mb)/p^ and M2 = {Mb - Ma)Ip^- 
Subtracting Eq. (|60|) from Eq. (|61|l similarly gives an equa- 
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tion for 5N, 



^ ^ on 



dt 



{AfiV^ + VA/i • V} 



6T 
d{5N) ■ 



(63) 



Equations H62|l and (|63|l can be cast into the more illu- 
minating form 



djSN) 



dt \ ^ Sn 



M-iS — 



A/iV 



(5((5iV)} 



(65) 



Equations (|64|l and H65|) couple the dynamics of the fields 
bN and n through a symmetric mobility tensor. The de- 
pendence of the mobilities Ma and Mb will in general 
depend on local crystal density and the local relative con- 
centration of species Ato B. 

For the case of substitutional diffusion between species 
A and B, Ma ~ Mb = M. In this limit the dynamics 
of n and bN decouple. Moreover if it is further assumed 
that that the mobility is a constant, Eas.f and H65f) 
become 



dn ^ ^2 
dt bn 



d{bN) 



AfeV^ 



bJ' 

b{bN) 



(66) 



(67) 



where the effective mobility Afg = 2M/ff. In the appli- 
cations using Eqs. and l|H7|l in the following sections, 
the dynamics of n and bN are simulated with time re- 
scaled hy t^t = 2Mt/p^. 

To illustrate the dynamics described by Eqs. (|66|l . H67|) 
and H43() . a simulation of heterogenous eutectic crystal- 
lization from a supercooled homogenous liquid was per- 
formed. The results of this simulation are shown in Fig. 
©. This figure shows the density (n), the density dif- 
ference {bN) and the local energy density at three time 
steps in the solidification process. These figures show 
liquid/crystal interfaces, grain boundaries, phase segre- 
gation, dislocations and multiple crystal orientations all 
in a single numerical simulation of the simple binary alloy 
PFC model. In this simulation a simple Euler algorithm 
was used for the time derivative and the spherical Lapla- 
cian approximation introduced was used (see Appendix 
A). The grid size was Aa; — 1.1 and the time step was 
Ai = 0.05. Unless otherwise specified all simulations to 
follow use the same algorithm, grid size and time step. 



V. APPLICATIONS 

This section applies the simplified phase field crystal 
model derived in Section IIII CI coupled to the dynam- 
ical equations of motion derived in Section IIVI to the 






^^^^^^^^^ 


^^^^^^ 


^^^^^^^^^^ 








FIG. 6: The grey scales in the figure correspond to den- 
sity (n), concentration (5N) and the local energy density in 
frames (a,b,c), (d,e,f) and (g,h,i) respectively. The area en- 
closed by white boxes is area shown in figures (a,b,c). The 
parameters in this simulation are the same as Fig. lO except 
L = 1.20 and Ri/Ro = 1/4 at JA^ = and ABo = 0.02. 
Figures (a,d,g), (b,e,h) and (c,f,i) correspond to times t = 
6600, 16200, 49900, respectively. 



study of elastic and plastic effects in phase transforma- 
tions. The first application demonstrates how the PFC 
alloy model can be used to simulate eutectic and den- 
dritic microstructures. That is followed by a discussion 
of the effects of compressive and tensile stresses in epitax- 
ial growth. Finally, simulations demonstrating the role 
of dislocations in Spinodal decomposition are presented. 



A. Eutectic and Dendritic Solidification 

One of the most important applications of the alloy 
phase field crystal model is the study of solidification 
microstructures. These play a prominent role in numer- 
ous applications such as commercial casting. Traditional 
phase field models of solidification are typically unable to 
self-consistently combine bulk elastic and plastic effects 
with phase transformation kinetics, multiple crystal ori- 
entations and surface tension anisotropy. While some of 
these effects have been included in previous approaches 
(e.g. surface tension anisotropy) they are usually intro- 
duced phenomenologically. In the PFC formalism, these 
features arise naturally from density functional theory. 

To illustrate solidification microstructure formation us- 
ing the PFC formalism, two simulations were conducted 
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FIG. 7: The grey-scale in the main portion of both figures 
show the concentration field 5N. In the insets, the grey-scale 
shows the density field, n, for the small portion of the main 
figure that is indicated by the white boxes, a) Eutectic crystal 
grown from a supercooled liquid at ABo = 0.0248 and 5N = 
0.0. The parameters that enter the model are the same as 
Fig. © except L = 1.20 and Ri/Ro = 1/4. b) Dendrite 
crystal grown from a supercooled liquid at ABo = 0.04 and 
SN — 0.0904. The parameters that enter the model are the 
same as Fig. ^ except L = 1.20 and Ri/Ro = 1/4. In Fig. 
(b) mirror boundary conditions were used. 



that the growth of a single crystal from a supercooled 
melt in two dimensions. In the first simulation a small 
perturbation in the density field was introduced into a 
supercooled liquid using the parameters corresponding 
to the phase diagram in Fig. jSJ, except L = 1.20 and 
Ri/Ro = 1/4. The reduced temperature ABq = 0.0248 
and average concentration of SN — 0.0. To reduce com- 
putational time the size of the lattice was gradually in- 
creased as the seed increased in size. A snapshot of the 
seed is shown at i = 480,000 in Fig. (7^). A simi- 
lar simulation was conducted for the growth of a den- 
drite from a supercooled melt for reduced temperature 
ABo = 0.04 and SN = 0.0904, with other parameters 
corresponding to those in Fig. ©, except for L — 1.20 
and Ri/Ro = 1/4. A sample dendritic structure is shown 
in Fig. Gt) at t = 175, 000. 

Simulations such as these can play an important role 



in establishing various constitutive relations for use in 
higher-scale finite element modeling (FEM) of elasto- 
plastic effects in alloys during deformation. In particular, 
traditional FEM approaches often employ empirical or 
experimental constitutive models to describe stress-strain 
response in elements that are intended to represent one 
(or more) grains. These constitutive relations are often 
limited in their usefulness as they do not self-consistently 
incorporate realistic information about microstructural 
properties that develope during solidification. 



B. Epitaxial Growth 

Another potential application of the PFC model is in 
the technologically important process of thin film growth. 
The case of heteroepitaxy, the growth of a crystalline film 
exhibiting atomic coherency with a crystalline substrate 
of differing lattice constant, has been examined in previ- 
ous PFC studies of pure systems 33, 34] . These initial 
works focused on two of the primary phenomena influenc- 
ing film quality: (i) morphological instability to buckling 
or roughening and (ii) dislocation nucleation at the film 
surface. A third important effect in alloy films, (iii) com- 
positional instability (phase separation in the growing 
film), requires consideration of multiple atomic species 
and their interaction. The purpose of this section is to 
illustrate how the binary PFC model addresses such com- 
positional effects in alloy heteroepitaxy, focusing on the 
spatial dynamics of phase separation over diffusive time 
scales. 

To date, a number of models of single component film 
growth incorporating surface roughening dislocation nu- 
cleation or both have been proposed j2^ |5^ [s^, [sl, 
Is^ l60l IbiL Ib^ l63j |. and models of binary film growth 
incorporating surface roughenin g an d phase separation 
have been proposed as well [63.1651 l6a . \6T\ . However, 
no existing models of binary film growth known to the 
authors have captured all of the above important phe- 
nomena, and it would be reasonable to expect that new 
insights into the nature of film growth could be gained 
through the simultaneous investigation of all of these 
growth characteristics. A unified treatment of this sort 
is required for the following reasons. There is clearly a 
strong link between surface roughening and dislocation 
nucleation, originating from the fact that dislocations 
nucleate at surface cusps when the film becomes suffi- 
ciently rough. It is also known that phase separation in 
the film is significantly influenced by local stresses, which 
are inherently coupled to surface morphology and dislo- 
cation nucleation. The dynamics of the growth process 
must then be influenced by the cooperative evolution of 
all three of these phenomena. In the next subsection 
numerical simulations will be presented to show that the 
binary PFC model produces all of the growth characteris- 
tics described above, and that each is influenced by misfit 
strain and atomic size and mobility differences between 
species. 
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1. Numerical simulations 

The physical problem recreated in these simulations is 
that of growth of a symmetric (50/50 mixture, SNq — 0) 
binary alloy film from a liquid phase or from a saturated 
vapor phase above the bulk coherent spinodal tempera- 
ture (ATc). Growth at temperatures above the miscibil- 
ity gap is typical of experimental conditions and should 
ensure that phase separation is driven by local stresses 
and is not due to spinodal decomposition. Initial con- 
ditions consisted of a binary, unstrained crystalline sub- 
strate, eight atoms in thickness, placed below a symmet- 
ric supercooled liquid of components A and B. In all 
the simulations presented, parameters are the same as in 
Fig. (01 except for L = 1.882 and AB = 0.00886 un- 
less specified in the figure caption. In what follows the 
misfit strain, e, is defined as [afum — asub)/asub, where 
o-fiim ~ o,a{^ + V^^o) if ill the constant concentration 
approximation. For a symmetric mixture of A and B 
atoms (i.e., 5No = 0) afUm (a^ + aB)/2. 

Periodic boundary conditions were used in the lateral 
directions, while a mirror boundary condition was ap- 
plied at the bottom of the substrate. A constant flux 
boundary condition was maintained along the top bound- 
ary, 120Aa; above the film surface, to simulate a finite 
deposition rate. Misfit strain was applied to the system 
by setting i? = 1 in the substrate and R = 1 -I- e -I- rjSN in 
the film. This approach yields a film and substrate that 
are essentially identical in nature except for this shift 
in lattice parameter in the film. Complexities resulting 
from differing material properties between the film and 
substrate are therefore eliminated, isolating the effects of 
misfit strain, solute strain, and mobility differences on 
the film growth morphology. The substrate was permit- 
ted to strain elastically, but was prevented from decom- 
posing compositionally except near the film/substrate in- 
terface. 

A sample simulation is shown in Fig. (jS)) for a mis- 
fit strain of e = 0.04, a solute expansion coefficient of 
r] = —1/4 and mobilities Ma = Mb = 1. As seen in 
this figure the well-documented buckling or Asaro-Tiller- 
Grinfeld j56l l57j instability is naturally reproduced by 
the PFC model. This instability is ultimately suppressed 
as a cusp-like surface morphology is approached, with 
increasingly greater stress developing in surface valleys. 
The buckling behavior ceases only when the local stress 
in a given valley imparts on the film a greater energy 
than that possessed by an equivalent film with a dislo- 
cation. At this stage a dislocation is nucleated in the 
surface valley and the film surface begins to approach a 
planar morphology. 

The nature of phase separation within the bulk film 
and at the film surface was found to vary with model pa- 
rameters, but a number of generalizations applicable to 
all systems studied have been identified. For the case 
of equal mobilities {Ma = Mb) we find that in the 
presence of misfit and solute strain, the component with 
greater misfit relative to the substrate preferentially seg- 




FIG. 8; Plots of the smoothed local free energy showing the 
progression of the buckling instability, dislocation nucleation 
and climb towards the film/substrate interface. From a) to f) 
times shown are t=150, 600, 1050, 1200, 1500 and 2550. In 
this figure e = 0.04, -q = -1/4 and Ma = Mb = 1. 



regales below surface peaks (see regions marked 1 and 
2 in Fig. ij^ and Fig. H10() '). Larger (smaller) atoms 
will be driven toward regions of tensile (compressive) 
stress which corresponds to peaks (valleys) in a compres- 
sively strained film and to valleys (peaks) in a film un- 
der tensile strain. This coupling creates a lateral phase 
separation on the length scale of the surface instabil- 
ity and has been predicted and verified for binary films 
[HlilililiiiiHEilzl and analogous be- 
havior has been predicted and verified in quantum dot 
structures fzlFlFsll. 

Secondly, again for the case of equal mobilities, the 
component with greater misfit relative to the substrate 
is driven toward the film surface (see Fig. Hl()|l ^. This 
behavior can also be explained in terms of stress relax- 
ation and is somewhat analogous to impurity rejection in 
directional solidification. The greater misfit component 
can be viewed as an impurity that the growing film wishes 
to drive out toward the interface. Experimental evidence 
from SiGe on Si ^] and InGaAs on InP 69, 70J verifies 
this behavior as an enrichment of the greater misfit com- 
ponent was detected at the film surface in both systems. 
Other models 0, EE 01 have not led to this type 
of vertical phase separation possibly due to neglecting 
diffusion in the bulk films. 

The third generalization that can be made is that, in 
the case of sufficiently unequal mobilities, the compo- 
nent with greater mobility accumulates at the film sur- 
face (see region marked 3 in Fig. ©). It was found that 
when the two components have a significant mobility dif- 
ference (typically greater than a 2:1 ratio) the effect of 
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FIG. 9: Plot of the smoothed local concentration field show- 
ing lateral phase separation between the surface peaks and 
valleys. White: Component A (large, fast), Black: Compo- 
nent B (small, slow). In this figure e = —0.02, rj = 0.4, 
Ma = 1, Mb = 1/4, and t = 3500. See text for discussion of 
the numbered arrows. 




FIG. 10: Plot of the smoothed local concentration field show- 
ing the nature of the phase separation under opposite signs of 
e. In figures a) and b) e = 0.04 and —0.025 respectively and 
in both figures Ma = Mb = 1 and rj = 0.25. 

mobility is more important than the combined effects of 
misfit and solute strains in determining which compo- 
nent accumulates at the surface. Since Ge is believed to 
be the more mobile component in the SiGe system, we 
see that the findings of Walther et al for SiGe on Si 
provide experimental support for this claim. They find 
a significant enrichment of Ge at the film surface, a re- 
sult that was likely due to a combination of this mobility 
driven effect as well as the misfit driven effect described 
in the second generalization. Experimental evidence also 
indicates that segregation of substrate constituents into 
the film may occur during film growth j?^ [tJ . We have 
similarly found that a vertical phase separation can be 
produced near the film/substrate interface and is com- 
plimented by a phase separation mirrored in direction 
near defects (see Fig. (|ll|l ). The extent of this phase 
separation is controlled largely by the bulk mobilities of 
the two constituents, and to a lesser degree by rj. The 
complimenting phase separation near climbing defects is 
a transient effect, any traces of which are dulled once the 
defect reaches the film/substrate interface. 



C. The Role of Dislocations in Spinodal 
Decomposition 

Spinodal decomposition is a non-equilibrium process 
in which a linearly unstable homogenous phase sponta- 
neously decomposes into two daughter phases. An ex- 
ample of this process in the solid state occurs during 



a) Tensile 









b) Compressive 



FIG. 11: Plot of the smoothed local concentration field show- 
ing the complimentary phase separation at the film/substrate 
interface and around defects. In this figure the film/liquid 
interface cannot be distinguished due to the overwhelming 
contrast created near the defects and at the film/substrate 
interface. In figures a) and b) e = —0.035 and 0.05 respec- 
tively, and in both figures Ma = Mb = 1 and rj = 1/4. 



a quench below the spinodal in Fig. Q when SN — 0. 
During this process domains of alternating concentration 
grow and coarsen to a tens of nanometers. Spinodal de- 
composition is of interest as it is a common mechanism 
for strengthening alloys, due to the large number of in- 
terfaces that act to impede dislocation motion. 

Solid state strengthening mechanisms, such as spinodal 
decomposition, rely critically on the interactions that ex- 
ist between dislocations and phase boundaries. Cahn was 
first to calculate that the driving force for nucleation of 
an incoherent second phase precipitate is higher on a dis- 
location than in the bulk solid 48]. A similar result 
was obtained by Dolins for a coherent precipitate with 
isotropic elastic properties in the solid solution '47'|. Hu 
et. al confirmed the results of Cahn and Dolins using 
a model that included elastic fields from compositional 
inhomogencities and structural defects |50l |. 

Recent studies of spinodal decomposition have used 
phase field models to examine the role of dislocations on 
alloy hardening [ioL l49l | . These phase field models couple 
the effects of static dislocations to the kinetics of phase 
separation. Leonard and Desai where the first to simu- 
late the effect of static dislocations on phase boundaries, 
showing that the presence of dislocations strongly favors 
the phase separation of alloy components [i^ . 

Haataja et al. recently introduced mobile dislocations 
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into a phase field model that couples two burgers vectors 
fields to solute diffusion and elastic strain relaxation. It 
was shown that mobile dislocations altered the early and 
intermediate time coarsening regime in spinodal decom- 
position 43, 78\ . Specifically, it was found that coherent 
strain at phase boundaries decrease the initial coarsening 
rate, since they increase stored elastic energy in the sys- 
tem. As dislocations migrate toward moving interfaces, 
they relax the excess the strain energy, thus increasing 
the coarsening rate l7a|. The growth regimes predicted 
by the model in Ref. |78j| are in general agreement with 
several experimental studies of deformation on spinodal 
age hardening [S El El Ei • 



1. Numerical Simulations 

The findings of Ref. ^7§\ were compared with numeri- 
cal simulations from the alloy PFC model. Here, spinodal 
decomposition was simulated for an alloy corresponding 
to the phase diagram in Fig. . Simulations began with 
a liquid phase of average dimensionless density difference 
SN — 0, which first solidified into a polycrystalline solid 
(alpha) phase, which subsequently phase separated as the 
reduced temperature (ABo) was lowered below the spin- 
odal. Figure (|12(l shows the concentration and density 
fields for four time sequences during the spinodal decom- 
position process. The dots in the figures denote the loca- 
tions of dislocation cores. Parameters for this simulation 
are given in the figure caption. Figure (|13() shows a plot 
of the average domain size versus time corresponding to 
the data in Figure fl^ . 

The PFC alloy results are in general agreement the re- 
sults of Ref. [73] ■ Specifically, PFC simulations show an 
early and intermediate time regime where the spinodal 
coarsening rate is is reduced from its traditional t^/"^ be- 
haviour. In this regime, the strain energy stored in the 
system is found to be much higher than that at late times 
when dislocations migrate to the boundaries, relaxing 
strain energy and leading to an an increased coarsening 
rate toward the usual t^^^ growth law. To more clearly 
illustrate the the interaction between a coherent bound- 
aries and dislocations, Fig. H14|) shows four time steps in 
the evolution of a dislocation migrating toward a phase 
boundary to relax mismatch strains. 

It is noteworthy that the alloy PFC introduced in this 
work does not incorporate "instantaneous" elastic relax- 
ation. A proper treatment of rapid relaxation of strain 
fields requires the model to be extended in a manner 
analogous to EJ. However, because of the asymptot- 
ically slow kinetics of spinodal decomposition and the 
small length scales between domain boundaries, it is ex- 
pected that this will only influence the time scales over 
which dislocations interact with domain boundaries. As 
a result, the general trends depicted in Figures. (|13|l and 
(|12|l are expected to be correct. 




FIG. 12: Four time sequences in the evolution of the concen- 
tration field (grey scale), superimposed on the corresponding 
density field. Dislocations are labeled by a square on the 
dislocation core surrounded by a circle. The time sequence 
(a)-(d) corresponds to t = 12000, 24000, 60000 and 288000, 
respectively (in units of At — 0.004). The system size is: 
1024Ax X 1024Aa;, where Ax = n/A. The density difference 
SN = 0, while L = 2.65, Ri/ = 0.25 and all other parame- 
ters are the same as Fig. 



VI. DISCUSSION AND CONCLUSIONS 

In this paper a connection between the density func- 
tional theory of freezing and phase field modeling was 
examined. More specifically it was shown that the phase 
field crystal model introduced in earlier publications 
[3^ 0, E3 and the regular solution commonly used in 
material science can be obtained from DFT in certain 
limits. These calculations relied on parameterizing the 
direct two-point correlation function that enters DFT by 
three quantities related to the elastic energy stored in 
the liquid and crystalline phases, as well as the lattice 
constant. 

In addition, a simplified binary alloy model was de- 
veloped that self-consistently incorporates many physi- 
cal features inaccessible in other phase field approaches. 
The simplified alloy PFC model was shown to be able 
to simultaneously model solidification, phase segrega- 
tion, grain growth, elastic and plastic deformations in 
anisotropic systems with multiple crystal orientations on 
diffusive time scales. 

It is expected that the alloy PFC formalism and its 
extensions can play an important role in linking material 
properties to microstructure development in a manner 
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FIG. 13: Inverse of the mean wave vector of the (circularly 
averaged) 2D structure factor of the concentration field Vs. 
time corresponding to the simulation in Figure 1121 

that fundamentally links the meso-scale to the atomic 
scale. As such this formalism, particularly when com- 
bined with novel adaptive mesh techniques in phase- 
amplitude or reciprocal space, can lead the way to a truly 
multi-scale methodology for predictive modeling of ma- 
terials performance. 
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FIG. 14: A dislocation migrates toward a coherent phase 
boundary thus relaxing mismatch strain. An 800 x 800 (units 
of Ax) portion of the actual simulations domain is shown. 
The data shows four time frames in the motion of the dis- 
location. Parameters of the simulation are the same as in 
Figure [H 

tions, t — uAt, X = iAx and y = jlS.x. 

An Euler discretization scheme was used for the time 
derivative and the 'spherical laplacian' approximation 
was used to calculate all Laplacians. For this method 
the discrete dynamics become 

V'ri+l.Jj = i'n,!,] + V^^„,jj, (Al) 

where pLn,i,j is the chemical potential given by 

^in^^,J = (r + (1 + V2)2)^„_,- ^. - , (A2) 
All Laplacians were evaluated as follows, 

+ fnS,i-l) 1'^ + (/n,i+lj + l + fnS-l,j + l 

+ j + 1 + j + l)/4 - 3/„,,j)/(Ax)2(A3) 



APPENDIX B: FREE ENERGY FUNCTIONAL 
FOR A PURE MATERIAL 



APPENDIX A: NUMERICAL DISCRETIZATION 
SCHEMES 

The governing Equations were numerically solved us- 
ing two different methods, described below. In what fol- 
lows the subscripts n, i and j are integers corresponding 
to the number of time steps and distance along the x and 
y directions, respectively, of a discrete numerical lattice. 
Time and space units are recovered by the simple rela- 



This appendix outlines the derivation of the free en- 
ergy functional T[p\ presented in Section 111 Al for a non- 
homogeneous liquid. The derivation follows the standard 
techniques of classical density functional theory as out- 
lined in 1541 (similar approaches are also employed in 
Ref. [1113^). 

The derivation begins with the definition of the grand 
partition function for N particles at temperature T, 

S = Tr{e(^"-''^)/^-«^} (Bl) 



17 



where 



Tr \ 



1 



(B2) 

is the classical trace operator, with fi and pt being the 
position momentum of the i*'' atom, respectively, fi is the 
chemical potential and h denotes Plank's constant. The 
N-body hamiltonian is denoted by = K + U + V, 
where 



N 



i=l 



N 



EL 



K 

U = U{n,--- ,nv) 
V 



(B3) 



i=l 



Here U denotes the interaction potential between par- 
ticles in the the system (including many body interac- 
tions), K is the total kinetic energy, with the mass 
of particles i and Vcxt represents the interaction of atom 
i with an external field. The probability density for any 
particular phase space configuration is given by 



/cq ^ 



(B4) 



The number density operator an N-body system is de- 
fined by by 



N 



(B5) 



The equilibrium number density is obtained by averag- 
ing the density operator with the equilibrium probability 
density. 



P{r) ^ (p) = Tr{/5/eq} 



(B6) 



The PFC formalism will ultimately yield governing equa- 
tions for the evolution of the number density in Eq. l|B6p 
on diffusive time scales. It is noteworthy that the equi- 
librium probability density, /eq, is a functional of p{f) 
|54| . an important property used in the derivation below. 
This is because for a given ?7, T4xt is uniquely determined 
by p{r)- Moreover, since Vcxt determines /oq, it follows 
that the equilibrium density /oq is a functional of p^r). 
The details of this argument are derived rigorously in in 
Ref. [s^l and will not be reproduced here. 

The arguments of the previous paragraph also imply 
that for a given [/, the Helmholtz free energy, defined by 



^[p] = Tr{/eq(i^ + 

is also a functional of p{f) 
tential functional, defined by 



f/ + fcsTln/eq)} (B7) 
I, as is also the grand po- 



n[p] = / drp{f)V,^tir) + ^[p] - p / drp{r) (B8) 



It should be noted that the grand potential can be cast 
into a more familiar form by substituting Eq. HB6() into 
f2[p] and exchanging the order of integration over r and 
Tr operation. Specifically, using the straightforward re- 
sults / drp{r)Vcxt{r) = -Tr{fcqV} and -p J drpif) = 
~Tr{fcqpN}, leads to the well known result from statis- 
tical mechanics, 



n[p] 



-ksT In ^ 



(B9) 



The grand potential is important in that it can be used to 
relate the chemical potential to the equilibrium density 
p{f) according to SV,[p] / 6 p{'r) = which gives [s^ 



P = Vcxt{r) + 



ST[p[ 
5p{r) 



(BIO) 



Equation (|B1Q() is fundamental to the theory of non- 
uniform fiuids in as it can, in principle, be used to calcu- 
late the equilibrium density |HlEl- 

The properties of the free energy functional T[p] can 
be be better elucidated by writing it as the sum of two 
terms 



(Bll) 



where Tq represents the ideal case of non-interacting par- 
ticles, while represents the total potential energy of 
interactions between the particles. Note that for a given 
[/, $ is, once again, a functional of p{r). Moreover, for 
[/ = in Eq. ||^ Tq[p\ becomes 



Mp] = keT / drp{r) (ln(Ap(f)) - 1) 



(B12) 



where A = {h? /2rmTkBT) 

To further discuss the properties of J'[p] for periodic 
phases, it will be useful to expand the free energy in 
Eq. (|B11(I about the density, p = pi, which corresponds 
to the liquid side of the solid-liquid coexistence phase 
diagram (at a given temperature). The change in free 
energy, Tc = J^[p] — J^[pi], then becomes 

= {To[p] - Toipi]) - mp] - ^pi]) (B13) 



The first term on the right hand side of Eq. (|B13|) can 
be simplified by substituting p ^ pi + Sp in the non- 
logarithmic expressions of Tq [p] , giving 

{To[p] - To[pi]) - keT f dr{p\n{p/pi) - 5p) (B14) 



In arriving at Eq. ljB14p . use was made of the property 
J dfSp = in the periodic state. The interaction term, 
(<i>[p] — '^[pi\)j can also be expanded functionally in 6p{r) 
about pi by defining 

5^p{r)] 



C2(ri,f2) 

C3(ri,r2,f3) 



6p{r) 



5p{ri)5p{f2) 
^p{'r3)Sp{ri)6p{f2) 



(B15) 
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Using Eq. ljB14|) and ljB15|l in Eq. ljB13|l finally gives values). Carrying out these expansions gives, 



dx 



p(f)ln(^)-5p(r) 



Pi 



-- J dridr25p{ri)C2{ri,r2)6p{r2) + 

-\ J dfidr2dr36p{fi)C3{ri,r2,r3)5p{r2)Sp{r3) 
6 

+ ••• (B16) 

The function C2 is the two point direct correlation func- 
tion of an isotropic fluid and it is usually denoted Cy = 
(^2(^12) , where ri2 = l^i — r2|. The function C3 is the 
three point correlation function, etc. 



J" 



pksT 



= / dr 



+ f (1 



-Caa + Cbb + '2Cab 



4 

,2 



Caa — Cbb 



({p-Pif + iP^-pD^ 



4p 

5N_ f _ Caa + Cbb - 2Cab 
2 V 4 



5N 



12 



where 



(C5) 



APPENDIX C: SIMPLE BINARY ALLOY MODEL 

This appendix goes through the expansion required to 
arrive at the simplified alloy model presented in section 
1111 Bl The starting point is the definition of the following 
three fields, 

n = {P-P)IP 
riA = ipA - Pa)/p 
riB = (pb - Pb)/p 

(CI) 

where overbars denote averages. The field n can equiva- 
lently be written as 

n = nA + nB ^ {pA + Pb)/p- {PA + Pb)/p (C2) 

Furthermore, the definitions of c = Pa/p and Sc = 1/2 — c 
(see Section lill B|) can be used to re- write 



/o = In i^-^j - (1 - Pi/p) ~ pC%/4 
-l{p + 2pll-p - ApdiC^A + C%b) 



(C6) 



Simplifying further, paying particular attention to the 
terms and substituting the explicit forms for C,"j gives, 



pkBT 



= / dr 



(3 



/o + B'y - y + ^ + nFV\ + nCV^n 



2p 

Caa — Cbb 
-pW)]5N 



+ 



SN 
~2~ 

6N^ 



{{P - Pif + {P" - Pl)n 



^ Caa + Cbb — '2Cab < 



-n + 



SN* 



12 



(C7) 



where 



^ PB - PA ^ pjriB - riA) + Pb- PA ^^^-^ 



2p 



2p 



which can further be used to define a new field 



SN = 2p5c/p ^ {UB - ua) + _ . (C4) 



Next, the fields p, Sp — p — pi, c and Sc in Eq. 
are expressed in terms of n and SN. Following that the 
free energy is expanded with respect to n and SN up to 
order four (noting that terms of order n or SN can be 
dropped since they integrate to zero in the free energy 
functional as they are all defined around their average 



= I- p{cu + c°BB + -^c^b)!^ 

C%u) 1 SN ^SN'' 



F = -p{C\a + CIb + 2Cls)/4 

^'^-iC\A-ClB)SN 
2p 

G = -p{C\a + C%b+2C\b)I^ 

+ f-{C\A-C%B)5N 
2p 



(C8) 



In what follows it is assumed that SN varies on length 
scales much larger than n. This is a reasonable on long 
time (diffusion) times scales, where solute and host atoms 
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intermix on length scales many times larger than the 
atomic radius. This assumption allows terms of order 
n to be eliminated from the free energy, i.e., 



where 



pksT 



fo 



+ 



df 

SN 
'~2~ 



2p 
Caa 



1 ?, A 

\ 

2 6 12 



1 



Caa + Cbb 



'2'Cab 



SN 



n^)SN 
Cbb 



nCW'^n 



12 



(^{i5-pif-pjn^))sN] (C9) 



and 



l-pCo 

^ + fjCo 
p 2p 



(C13) 



The previous equation can finally be cast into a form 
similar to that presented in Section IlII CI of the text, 



= / dr 



S^(2i?2v2+i?4v4)j 



n ui 
12 ^ 2" 



SN^ 
IT 



|V(5iV|^ 



+ -/SN + ^SNV^SN 



(CIO) 



where 



R 

w 
L' 

7 



FV(2G) 
V2GJf 
(1 - {C 

{C\ 
~{C 



'AA + 
^AA 



'0 

AA 
'^BB 



c° 



2CU)/2) 



^BB 



BB 

^ciB)m 

2C1b)/2) 



2p 

Caa 



(l-r.3) 



Cbb ( 



A-p 



pW 



(Cll) 



The dependence of the coefficients in _B', B^ and R 
on the density difference can be explicitly obtained by 
expanding them in SN as well. This gives, 



B' = 

= 

R = Rn 



B[ SN + B^ SN'^ 
BISN + B^ SN' + ■■■ 
Ri SN + R2SN' + ■■■ 



(C12) 



while 



and 



Cn 

SC., 



(c 



AA 



^ AA ^BB-> 



C^bI + ^C^aI)/^ 

(n) 



Bl = 
BH = 



SC,f 
Ca 



c 



IP I 



ApCi 



{C2SC4 — 2CiSC2) 



-XC^SCA-C^SC^f 



Ro — 
Ri/Rq — 
R2/R0 = 



2C4 
\ 4 



pI 



(C14) 



(C15) 



Ap'CiC2 



-{C2SC4-C4SC2) 



32piCjCi 

{d2SCi + ?,CiSC2) 



((72(5(74 - (74(5(72) X 
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